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Abstract 

An axisymetric phase field model is developed and 
used to model surface tension forces on liquid jets in 
microgravity. The previous work in this area is 
reviewed and a baseline drop tower experiment selected 
for model comparison. A mathematical model is 
developed which includes a free surface, a symmetric 
centerline and wall boundaries with given contact 
angles. The model is solved numerically with a 
compact fourth order stencil on a equally spaced 
axisymetric grid. After grid convergence studies, a grid 
is selected and all drop tower tests modeled. Agreement 
was assessed by comparing predicted and measured 
free surface rise. Trend wise agreement is good but 
agreement in magnitude is only fair. Suspected sources 
of disagreement are suspected to be lack of a turbulence 
model and the existence of slosh baffles in the 
experiment which were not included in the model. 

Nomenclature 

We = Weber Number 

p = density 

u = velocity 

r =jet radius 

a = surface tension 

f = free energy 

a = constant 1 

(3 = constant 2 

C = phase distribution 

v|; ^barrier function 

V = volume 

(j) = potention function 

P = pressure 

x = distance 

g = wall function 

M = radial correction factor 

9 C = contact angle 

a ^height 

d = jet diameter 

D = tank dimater 

L = tank length 

Re = Renolds Number 


Introduction 

Microgravity poses many challenges to the designer of 
spacecraft tanks. Chief among these are the lack of 
phase separation in the fluid and the need to supply 
vapor-free liquid or liquid free vapor to the required 
spacecraft processes. One of the principal problems of 
phase separation is the creation of liquid jets. A jet can 
be created by liquid filling, settling of the fluid to one 
end of the tank, or even closing a valve to stop the 
liquid outflow. In normal gravity the gravitational force 
controls and restricts the liquid jet flow, but in 
microgravity with gravity largely absent, jets must be 
contained by surface tension forces. Recent NASA 
experiments in microgravity (TPCE and VTRE) have 
brought a wealth of data of jet behavior in 
microgravity. VTRE was surprising in that although it 
contained a complex geometry of baffles and vanes the 
limit on liquid inflow was the emergence of a liquid jet 
from the top of the vane structure. Clearly 
understanding the restraint of liquid jets by surface 
tension is key to managing fluids in low gravity. 

Flow of a submerged axial jet constrained by surface 
tension is similar to stagnation flow against a plate in 
that the jet hits the constraining surface and is deflected 
radially out. However, the ability of the constraining 
surface to move in response to the exerted force is 
unique. In fact to increase the restraining force on the 
jet as flow rate increases, the surface must deform to 
increase the radius of curvature of the free surface, 
thereby increasing the surface tension force. 
Unfortunately this motion also changes the direction in 
which the surface tension force acts. Eventually the 
limit is reached where the surface tension force is 
perpendicular to flow and hence no longer retains it. 
When the deformation of the free surface is large the 
restraining bulge is long and slender. At this point 
several other mechanisms act to break down the jet, 
such as columnar buckling or the Taylor Instability 
where surface waves grow to such amplitude that they 
pinch a droplet off from the jet. 
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To model this phenomenon a numerical method that 
tracks the fluid motion and the surface tension forces is 
required. Jacqmin 1 has developed a phase-field model 
that converts the delta-function surface tension force 
into a continuum function that peaks at the free surface 
and decays rapidly away. Previous attempts at this 
formulation have been criticized for smearing the 
interface but by sharpening the phase function, double 
gridding the fluid function and using a higher order 
solution for the fluid function these concerns have been 
ameliorated. 

Review of Literature 

NASA Drop Tower Data is found in references 2-13. 
Symons and Spuckler 7 studied the liquid inflow via 
axial jet into a broad range of tank shapes both empty 
and partially full. Symons work establishes an 
empirical limit for jets of Weber number (We) equal to 
l.j-1.5 depending on jet velocity profile, where: 


p — density 

u av = average jet velocity 
r = jet radius 
a = surface tension 

Staskus extends the work of Symons by placing 
baffles in front of the jet. However, no attempt is made 
to analyze these complex flows. Instead results are 
reported as a ratio of improvement to the unbaffled jet 
Weber number. Labus 8 also studies the effect of baffles 
including ones that break the central jet into several 
small jets. Aydelott 10 ’ 12 ' 13 looks at the problem of a 
recirculating jet where the liquid level is held constant. 
Results are classified into four flow patterns, 
dissipation, geyser formation, aft collection, and 
circulation. It is the geyser formation/aft collection we 
concern ourselves with in this paper. Aydelotfs 
assessment that a drop in mixing accompanies this 
transition indicates the transition's importance. Labus 11 
studies both stagnation flow and free surface shape, but 
is concerned with the back of a flow stagnated against a 
flat plate in microgravity. 

Shuttle based experiments in references 14-18 provide 
valuable flight data. Video of Plexiglas tanks during 
shuttle flight provide several improvements over drop 
tower tests; including increasing the scale from 4” tanks 
to 12” tanks and extending the duration of test from 5 
seconds to half-an-hour. Tank Pressure Control 


Experiment is the first of these and has flown three 
times. The first flight focused on the mixing studies of 
Aydelott. Improvements included actual heat transfer 
data by using a condensing fluid (refrigerant 1 13) and 
longer duration. Bentz 14-16 was able to confirm the 
geysering and circulating regimes of Aydelott, but 
encountered an asymmetric regime between the two 
that was even more catastrophic to heat transfer than aft 
collection. The second flight of TPCE focused mostly 
on rapid boiling phenomena, but contains some further 
tests on mixing. Hasan 17 confirms the findings of 
Bentz. The third flight 18 was done at a lower fill level 
but confirms the results of the other flights. The Vented 
Tank Resupply Experment 19 was designed to look at 
vanes rather than axial jets, but as noted previously 
exhibits the classic gesyering behavior. 

Analytical work is listed in references 20-30. Concus 20 ' 
provides differential equations of the free surface 
problem, but analyzes only static cases. Nickell 22 
analyzes flow from a jet into a liquid and the resultant 
free surface shape for a normal gravity application, but 
removes all surface tension from the analysis as 
secondary. Hochstein 23 ' 24 analyzes the microgravity 
mixing with a volume of fluid approach, but uses only a 
limited approximation to model the surface tension. 
Aydelott 25 and Der 26 both analyze the motion of a 
bubble in the oxygen tank during separation of a 
Centaur stage with VOF models; noteworthy in these is 
again the appearance of a geyser. Tegart 27 shows the 
application of the surface Evolver code of Brakke 28 to 
actual tank shapes. BrackbilP 9 develops an improved 
surface tension model for VOF codes, but only shows 
one example of its use for axial jets. Shrader 30 uses a 
Runga-Kutta scheme to solve the differential equation 
of free surface deformation in response to an imposed 
pressure field. This approach is quite promising but 
does not always converge and limits the interaction 
between the flow field and the free surface. Jacqmin 1 
developed a phase field model of surface tension and 
implemented as a fourth order accurate scheme using a 
compact 9-point stencil. Although Jacqmin lays out the 
basic axisymmetric scheme the computer code and all 
the examples in his paper are planar. The Jacqmin 
model will serve as the basis for the present analysis 

Model 

Introduction 

To model the fluid motion the Navier-Stokes equations 
are formulated for low-speed incompressible flow. 
Velocity and pressure are placed on a staggered grid, 
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with velocity being tracked at cell faces and pressure at 
cell centers. To track the free surface a color function is 
introduced which tracks liquid as 1/2 and gas as -1/2. 
Jacqmin has developed a phase-field model that 
converts the delta-function surface tension force into a 
continuum function that peaks at the free surface and 
decays rapidly away. Previous attempts at this 
formulation have been criticized for smearing the 
interface but by sharpening the phase function, double 
gridding the fluid function and using a higher order 
solution for the fluid function these concerns have been 
ameliorated. This paper will document the adaptation 
of the Jacqmin algorithm to the problem of restraint of 
liquid jets. The enhancements include formulation of an 
axisymmetric fourth order model, implementation of a 
symmetric boundary condition at the tank centerline, 
and extension of the wall wetting boundary condition to 
fourth order accuracy. A simple velocity forcing 
function has been added to simulate the jet without 
violating continuity. 

Phase Model of Surface Tension 

Surface tension can be expressed as a free energy field. 
The expression for this energy in our formulation is 
given by 

f = \a\VC'\ +(W{C) (2) 

Where C is a phase distribution function and vp is a 
barrier function that is maximum at the interface and 
dies away as the phase becomes uniform. This 
formulation is extracted from Van der Waals 31 and 
inherently impliesthat the equilibrium free surface 
position is the one where the free energy is minimized. 

In order to model this behavior the physical vp that dies 
away on the molecular scale is approximated by a 
function with similar behavior on a larger scale such as 


,2* 


y(C) = 


2k + 2 


■C 


2k+2 


- C 2 + — 

2 S(2k + 2) 


(3) 


This function has the required properties of being 
maximum at C=0 and dying away to 0 at both 1/2 and - 
1/2. If we define our C function as being 1/2 when the 
phase is liquid and -1/2 when the phase is gas this will 
produce the required behavior. Higher values of k 
produce sharper peaks. For our solution we will choose 
k = 16. 


To study the transients of the free surface some 
additional formulations are required. We define a 
potential function as the rate of change in f per unit 
volume with respect to C. 


5 \fdV 

0 = = /W'(C)-aV 2 C (4) 

o C 

Cahn and Hilliard j2 approximate the transients of the 
free surface by setting the diffusion fluxes as 
proportional to the potential gradient. In equation form 
this is 

f = ( 5 ) 

This gives us two coupled Possion equations to solve 
for the phase distribution. To add the effects of fluid 
motion we must use the Navier-Stokes equations. The 
continuity equation for incompressible flow is 


V • u = 0 (6) 

The momentum equations for each direction are given 
by 


Du, du, sr- du. 

P^7 = P-^ + r > T. u iir= 


Dt ' dt 

72. 


j fej (7) 
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Fourth Order Formulation of the Governing 
Equations 


The equations of the previous section cannot be solved 
directly but must be solved numerically. To keep the 
interface as sharp as possible a compact 4 th order stencil 
is used 

The mehrstellungen method 
The equation 

V 2 u = f (8) 

Was show by Collatz 33 to be approximated by 

W + 1(a*) ! v‘'U/+.1v ! / 

L 12 J 12 (9) 

+0( Ax 4 ) 


Using central differencing on a square Cartesian grid 
one obtains the following computational stencil 
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Wall Boundaries 


With some slight modification we can use this to 
rewrite the potential equation as 

n 4 ii 
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We retain the V 2 <j> term because it can be calculated by 
equation 5. 

Application to an Axisymmetric Grid 
In axisymmetric coordinates our equation becomes 
’ M_ 4 M + ~ 

4 M_ -20 4M C = 
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where 


M ± = 0 ± +r 0 )/2r 0 


Symmetry Boundary 

At the center line (r=0) the problem is symmetric 
Theretore along this edge the equation becomes: 
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At the outer wall two boundaiy conditions will be used. 
First is the no flux boundary 

30 n 

¥ =0 <15) 
The second is a bit more complicated. Postulating a 
wall energy function of the form 

fw = j&g(C)dA (16) 

where g is a function chosen to yield the correct contact 
angle, then the diffusively controlled equilibrium at the 
wall is 

a— + c7g\C) = 0 (17) 

dr 

for our purposes 

g(C) = cos(# c )[6C 2 -1.5] ( 18 ) 

where 0 C is the static contact angle at the side wall. 

For the experiment modeled the contact angle is zero so 
COS(f? ) is one. 

Substituting into our main equation 
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A similar approach is used for the top wall except that 
for simplicity 0 C is set to 90° resulting in g=g — 0 and 
the boundary being symmetric. This should not affect 
the results of the calculation since interface flow along 
the top wall does not occur in our problem until after 
free surface penetration. The equation for the top wall 
is 
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The bottom wall is done the same as the top wall. Here 
the logic used to justify 0 C is set to 90° is that the wall is 
only in contact with liquid throughout our runs. The 
equation for the bottom wall is 
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Comer Boundaries 

Equation for the top inner boundary combining the 
symmetry and top wall boundaries 
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Equation for the bottom inner boundary combine 
symmetry and bottom wall boundary 
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Equation for the top outer boundary combining top and 
outer boundary 
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Equation for the bottom outer boundary 
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Implementation as a CFD code 

The previous equations form a complete set of 
differential equations that can be solved for the fluid 
transient motion. Each equation is solved sequentially 
and numeric techniques specific to each equation are 
used to achieve the desired level of accuracy. 

Solution for the potential field 

A Newton-Rhapson iteration is used to project the 
body centered values of C to the cell boundaries and 
produce a c|> field consistent with equation (5) 

Advection of phase quantities 

The equation (12) and its boundary equations form a 
matrix equation that is solved using the current 
values of C and (j> to project new values. This process 
is iterated four times to smooth the solution. 


Solution of the velocity field 

Equation (7) is used to predict the change in velocity 
field. The projected velocity changes are used to 
calculate viscous stresses that are then used to correct 
the velocity change. 
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Solution of the pressure equation 

The velocity changes are fed into the pressure Poisson 
equation that is solved by successive over relaxation to 
produce a uniform static pressure field consistent with 
our incompressible flow assumption. 

Approximation of the Liquid Jet 

Since the solution of the Navier-Stokes equation we 
used conserves mass strongly creating source and 
sink flows pose many difficulties. To avoid these 
problems the entering jet was modeled as a 
recirculating region where the v velocity was forced 
to a desired value. This allowed the u velocity to 
entrain liquid into the jet and thus conserve mass. A 
0.5 cm length for this region was chosen since this 
appeared to be long enough so at the top of the region 
the mass entrainment was sufficient such that the v 
velocity was the dominant fluid motion. 


Comparison to Experiment 
Model Runs of Test Cases 

After implementation of the code in axisymmetric form 
was complete and verified by several test cases, the 
drop tower runs of Aydelott were modeled. Aydelott 
looked at the problem of a recirculating jet where the 
liquid level is held constant. Results are classified into 
four flow patterns, dissipation, geyser formation, aft 
collection, and circulation. It is the dissipation/geyser 
formation we concern ourselves with in this paper. 
Table 1 shows a compilation of Aydelott’ s zero-g runs 
in these regimes. Little of the drop tower film remains, 
but figure 1 shows the time history of run 15 (Re 450, 
Fill level 50%). Times are estimated from frame counts 
since the clock is out of focus. Spherical tank data has 
been omitted since we cannot model curved 
boundaries. Aydelott' s assessment that this transition is 
accompanied by a drop in mixing indicates the 
transition’s importance. As a turbulence model for this 
formulation is not yet available the work of Aydelott 
was deemed the most suitable since it contains a 
number of cases observed by the author to be either 
laminar or transitional with Reynolds number between 
450 and 1290. Figures 2-4 show representative 
computer predictions of these tests. Figure 2 is a low 
flow rate test that only slightly deforms the free surface. 
Figure j is a medium flow rate test that forms a geyser 
in the center but reaches a steady flow condition. Figure 
4 is a high flow rate test where the geyser continues to 
grow throughout until it eventually contacts the far 
boundary of the grid. The model handed the free 


Table 1 Experimental Results of Aydelott 


Test 

Tank 

Shape 

Liquid 
Fill 
Vol % 

Jet 

Reynolds 

no 

Jet 

Weber 

no 

Ratio of 
Geyser 
height to 
tank 

Diameter 

1 

c 

29 

630 

0.96 

0.55 

2 

c 

29 

900 

1.09 

.80 

12 

b 

39 

450 

.49 

.36 

13 

b 

39 

900 

1.11 

.84 

14 

b 

39 

1290 

1.59 

2.16 

15 

c 

50 

450 

.39 

.42 

17 

c 

51 

630 

.78 

.34 

19 

c 

52 

900 

.81 

.52 

24 

a 

52 

1320 

1.16 

1.45 

50 

b 

60 

450 

.37 

.24 

51 

b 

60 

900 

.72 

.42 

52 

b 

60 

1320 

1.05 

1.10 

53 

b 

73 

900 

.57 

.30 

57 

b 

73 

1270 

.78 

.70 

64 

b 

91 

480 

.31 

.10 

65 

c 

91 

900 

.48 

.20 

66 

c 

92 

1290 

.62 

.48 


surface deformation quite nicely, even to the point of 
modeling geyser growth in the regime where the free 
surface is no longer restrained. 


Grid sensitivity 

A grid sensitivity test was done to confirm the choice of 
a 50 by 200 grid for modeling. This grid is fine enough 
to place two points in the starting jet, but yet not overly 
tax the computer for storage and run time. Comparison 
to a 100 by 400 grid showed no change in either 
flowfield or free surface shape. Figure 5 shows a 
double grid run comparable to figure 2. 

On Laminarity of the Jet Inflow 

Classic analysis of the liquid in liquid jet indicates very 
low stability. Reference 35 shows a limit as low as Re 
= 11. Mcnaughtoun and Sinclair using a more practical 
analysis divided the liquid in liquid jet into four 
regimes 

dissipated laminar (Re < 300 approx.) 

Fully laminar jets (300 < Re < 1000 approx.) 
semi-turbulent (1000 < Re < 3000 approx) 

Fully turbulent (Re > 3000) 

Their transition numbers are somewhat a function of 
the length of their apparatus. Dissipated jets only made 
it a portion of the way across the test chamber before 
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dissipating into the bulk liquid via the same breakdown 
reported in the previous reference. This distance 
increased with increasing Re number until a laminar jet 
spanned the length of their test chamber the fully 
laminar region. As the Re increased further a turbulent 
flow region began to emerge near the far end of the jet. 
The length of laminar flow would begin to decrease as 
the flow increased until eventually the jet would 
become turbulent right at the nozzle marking the 
transition to the final region. 

Their data for the length of the laminar region was 
given by the correlation 

ald = 9.97x1 0 7 Re' 2 46 (£> / d)~°" 48 ( L / df 74 (26) 

Our runs of Re 450, 630 are in the fully laminar region. 
The Re 900 run should transition to turbulence at a 
distance of 8.3 cm and Re 1290 at 3.4 cm. Aydeiott 
reviewing the data reports the following findings. No 
spreading at Re 450, for Re 630 no spreading if the 
liquid height over the jet was less than 2.5 cm; 
spreading consistent with a laminar jet thereafter. Re 
900 a jet intermediate between laminar and turbulent. 

And for Re > 1500 a fully turbulent jet. Comparison to 
drop tower data 

Visual Comparison 

Comparison of the data to the model show similarity in 
jet spread and flow motion. The model even captures 
the vortex shedding from the tip of the geyser as the 
flow develops although the axisymetric nature of the 
model forces more regularity in the vortex shedding 
than is seen in the drop tower film. 

Predicted Geyser Height 

Model predictions of geyser heights are shown in table 
2. For comparison the measured heights of Aydeiott 
and the previous predictions of Schrader 30 (Only 
available for a few runs) are shown. 

Figure 6 shows the comparison between predicted 
geyser height and fill level for the model and the 
experiment. Although the model underpredicts 
experiment for the lower flow rates it overpredicts at 
the highest flow rate. The overprection is believed due 
to the lack of turbulence modeling. It is well known 
that turbulent jets spread at a much higher rate than 
laminar jets. This increased spread will lower the 
centerline velocity more quickly and increase the area 
of the jet at the free surface, decreasing the amount of 
surface deformation required to contain the jet. The 
under prediction is harder to explain but is worse at low 
fill levels. It may be due to a slosh baffle at the 33% full 
level which prevents development of a hemispherical 


Table 2 Geyser Height Comparison 


Test 

Model 
Ratio of 
Geyser 
height to 
tank radius 

Prediction 
of Schrader 

Measured Ratio 
of Geyser height 
to tank radius 

1 

0.22 

N/A 

0.55 

2 

0.4 

N/A 

.80 

12 

0.12 

N/A 

.36 

13 

0.42 

0.36 

.84 

14 

2.7 

N/A 

2.16 

15 

0.06 

N/A 

.42 

17 

0.1 

N/A 

.34 

19 

0.28 

N/A 

[“52 

24 

2.06 

N/A 

1.45 

50 

0.04 

N/A 

.24 

51 

0.1 

N/A 

.42 

52 

1.84 

N/A 

1.10 

53 

0.2 

0.21 

.30 

57 

0.84 

0.45 

.70 

64 

0 

N/A 

.10 

65 

0.16 

N/A 

.20 

66 

0.78 

N/A 

.48 


interface in zero gravity. The baffle acts to raise the 
liquid height at the centerline and flatten the free 


surface which lessens the surface tension force on the 
jet. This effect was most pronounced at the 29% full 
runs where the measured height of the free surface 
above the jet is actually higher than the 39% full tests. 
It was felt that free suface height above the jet was the 
more important parameter. Therefore model runs 


matched the free height not the liquid fill level. 


Overall agreement is reasonable. Run 14 will be 
difficult to match since this run is most likely unstable 
so the measured geyser height is limited by the tank 
size. Run 13 is as yet unexplained but Schrader’s 
prediction is even more in error. For the other two runs 
predicted by Schrader, comparison is as follows. Run 
53 both predictions are very close to each other but 
low. Run 57 the current prediction is much closer to the 
measured value. 


Summary 

This paper documents the adaptation of the Jacqmin 
algorithm to the problem of restraint of liquid jets. 
Adaptations include formulation of an axisymmetric 
fourth order model, implementation of a symmetric 
boundary condition at the tank centerline, and 
enhancement of the wall wetting boundary condition to 
fourth order accurate. A simple velocity forcing 
function has been added to simulate the jet without 
violating continuity. Figures demonstrate the code’s 
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ability to model jet flows. Comparison to drop tower 
data shows the strength and validity of this code. 
Finally, the limits of surface tension in restraining 
liquid jets have been found and the growth rate of 
geysers after this limit has been exceeded calculated. 
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